###############################################################################   
#### Replication Materials                                                 #### 
#### Kim, Nakka, Gopal, Desmrais, Mancinelli, Harden, Ko, Boehmke. 2021.   ####
#### Attention to the COVID-19 pandemic on Twitter:                        ####
#### Partisan differences among U.S. state legislators                     ####
#### Legislative Studies Quarterly                                         ####
###############################################################################  


###############################################################################
################################### Set Up ####################################
###############################################################################

# packages -------------------------

lapply(c('readr', 'stargazer', 'texreg', 'xtable', 'lme4', 'glmmTMB','plm',
  'lmtest', 'interplot','dummies','effects', 'ggthemes','caret', 
  'ggeffects', 'scales','ggrepel','sandwich','reshape2'), 
  require, 
  character.only = TRUE)

# read data set for regression -------------------------

data_path <- '/Users/taegyoon/Google Drive/spap_state/spap_state_attention/data/' 
df <- read_csv(paste0(data_path, "spap_state_attention_regression.csv")) 
df <- subset(df, select = -c(X1))


###############################################################################
################## Main Regression Models: Table 1, Table S5 ##################
###############################################################################

# set the panel structure -------------------------

plm_df <- pdata.frame(x = df, 
                      index = c('user.screen_name', 'week'))

# fit models -------------------------

main_state <- plm(covid_relevant_1_log ~ 
                    state_case_pop + 
                    state_death_pop + 
                    state_covid_policy + 
                    party_new + 
                    state_abbrev + 
                    week_new + week_new_2 + week_new_3,
                  data = plm_df,
                  index = c('user.screen_name', 'week_new'),
                  model = 'random', 
                  effect = 'twoways')

main_national <- plm(covid_relevant_1_log ~ 
                       national_case_pop + 
                       national_death_pop + 
                       state_covid_policy + 
                       party_new + 
                       state_abbrev + 
                       week_new + week_new_2 + week_new_3,                                             
                     data = plm_df,
                     index = c('user.screen_name', 'week_new'),
                     model = 'random', 
                     effect = 'twoways')

main_all <- plm(covid_relevant_1_log ~ 
                  state_case_pop + 
                  state_death_pop +
                  national_case_pop + 
                  national_death_pop +
                  state_covid_policy + 
                  party_new +
                  state_abbrev +
                  week_new + week_new_2 + week_new_3,            
                data = plm_df,
                index = c('user.screen_name', 'week_new'),
                model = 'random', 
                effect = 'twoways')

main_state_inter <- plm(covid_relevant_1_log ~ 
                          party_new*state_case_pop + 
                          party_new*state_death_pop +
                          party_new*state_covid_policy + 
                          state_abbrev +
                          week_new + week_new_2 + week_new_3,            
                        data = plm_df,
                        index = c('user.screen_name', 'week_new'),
                        model = 'random', 
                        effect = 'twoways')

main_national_inter <- plm(covid_relevant_1_log ~ 
                             party_new*national_case_pop + 
                             party_new*national_death_pop +
                             party_new*state_covid_policy + 
                             state_abbrev +
                             week_new + week_new_2 + week_new_3,            
                           data = plm_df,
                           index = c('user.screen_name', 'week_new'),
                           model = 'random', 
                           effect = 'twoways')

main_all_inter <- plm(covid_relevant_1_log ~ 
                        party_new*state_case_pop + 
                        party_new*state_death_pop +
                        party_new*national_case_pop + 
                        party_new*national_death_pop +
                        party_new*state_covid_policy + 
                        state_abbrev +
                        week_new + week_new_2 + week_new_3,            
                      data = plm_df,
                      index = c('user.screen_name', 'week_new'),
                      model = 'random', 
                      effect = 'twoways')

main_all_inter_ctl <- plm(covid_relevant_1_log ~ 
                            party_new*state_case_pop + 
                            party_new*state_death_pop +
                            party_new*national_case_pop + 
                            party_new*national_death_pop +
                            party_new*state_covid_policy + 
                            np_score +
                            major_cham +
                            state_abbrev +
                            week_new + week_new_2 + week_new_3,            
                          data = plm_df,
                          index = c('user.screen_name', 'week_new'),
                          model = 'random', 
                          effect = 'twoways')

# adjust SE for serial correlation -------------------------

main_state_adj <- coeftest(main_state, 
                           vcovHC(main_state, 
                                  method = 'arellano'))

main_national_adj <- coeftest(main_national, 
                              vcovHC(main_national, 
                                     method = 'arellano'))

main_all_adj <- coeftest(main_all, 
                         vcovHC(main_all, 
                                method = 'arellano'))

main_state_inter_adj <- coeftest(main_state_inter, 
                                 vcovHC(main_state_inter, 
                                        method = 'arellano'))

main_national_inter_adj <- coeftest(main_national_inter, 
                                    vcovHC(main_national_inter, 
                                           method = 'arellano'))

main_all_inter_adj <- coeftest(main_all_inter, 
                               vcovHC(main_all_inter, 
                                      method = 'arellano'))

main_all_inter_ctl_adj <- coeftest(main_all_inter_ctl, 
                                   vcovHC(main_all_inter_ctl, 
                                          method = 'arellano'))


###############################################################################
################### Generate Latex Code : Table 1, Table S5 ###################
###############################################################################

# get information for the last six rows -------------------------

texreg(
  list(main_state, main_national, main_all,
       main_state_inter, main_national_inter, main_all_inter, 
       main_all_inter_ctl),
  omit.coef = 'abbrev',
  digits = 5,
  caption = 'Panel Regression Model (State Fixed Effect and Legislator-week Random Effect): Population Normalized + Three-fold Party Variable + Week Polynomial + Logged DV (1 added)',
  custom.coef.names = c('(Intercept)',
                        'State New Cases (per 10k)',
                        'State New Deaths (per 10k)',
                        'State COVID-19 Policies',
                        'Other',
                        'Republican',
                        'Week',
                        'Week (quadratic)',
                        'Week (cubic)',
                        'National New Cases (per 10k)',
                        'National New Deaths (per 10k)',
                        'Other * State New Cases (per 10k)',
                        'Republican * State New Cases (per 10k)',
                        'Other * State New Deaths (per 10k)',
                        'Republican * State New Deaths (per 10k)',
                        'Other * State COVID-19 Policies',
                        'Republican * State COVID-19 Policies',
                        'Other * National New Cases (per 10k)',
                        'Republican * National New Cases (per 10k)',
                        'Other * National New Deaths (per 10k)',
                        'Republican * National New Deaths (per 10k)',
                        'Legislator Ideology',
                        'Chamber Majority Status'
                        )
  )

# get information for betas and corresponding standard errors -------------------------

texreg(
  list(main_state_adj, main_national_adj, main_all_adj,
       main_state_inter_adj, main_national_inter_adj, main_all_inter_adj, 
       main_all_inter_ctl_adj),
  omit.coef = 'abbrev',
  digits = 5,
  caption = 'Panel Regression Model (State Fixed Effect and Legislator-week Random Effect): Population Normalized + Three-fold Party Variable + Week Polynomial + Logged DV (1 added)',
  custom.coef.names = c('(Intercept)',
                        'State New Cases (per 10k)',
                        'State New Deaths (per 10k)',
                        'State COVID-19 Policies',
                        'Other',
                        'Republican',
                        'Week',
                        'Week (quadratic)',
                        'Week (cubic)',
                        'National New Cases (per 10k)',
                        'National New Deaths (per 10k)',
                        'Other * State New Cases (per 10k)',
                        'Republican * State New Cases (per 10k)',
                        'Other * State New Deaths (per 10k)',
                        'Republican * State New Deaths (per 10k)',
                        'Other * State COVID-19 Policies',
                        'Republican * State COVID-19 Policies',
                        'Other * National New Cases (per 10k)',
                        'Republican * National New Cases (per 10k)',
                        'Other * National New Deaths (per 10k)',
                        'Republican * National New Deaths (per 10k)',
                        'Legislator Ideology',
                        'Chamber Majority Status'
                        )
  )